Density functional theory with adaptive pair density 
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We propose a density functional to find the ground state energy and density of interacting particles, 
where both the density and the pair density can adjust in the presence of an inhomogeneous potential. 
As a proof of principle we formulate an a priori exact functional for the inhomogeneous Hubbard 
model. The functional has the same form as the Gutzwiller approximation but with an unknown 
kinetic energy reduction factor. An approximation to the functional based on the exact solution of 
the uniform problem leads to a substantial improvement over the local density approximation. 
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Most of our theoretical understanding of condensed 
matter and complex molecules stems from density func- 
tional theory (DFT) computations [1]. Practical im- 
plementations rely on the local density approximation 
(LDA) or its gradient generalizations which often pro- 
vide accurate results at the modest cost of a Hartree 
like computation [2| . These methods however fail in sys- 
tems at or close to the Mott insulating regime [3[, when 
the tunneling matrix element of electrons becomes small 
compared with the typical electron-electron repulsion en- 
ergies. This can be seen already at the level of an H2 
molecule which is stretched to produce two separated H 
atoms. LDA performs reasonably well at the equilib- 
rium distance, when electrons have substantial tunneling 
among the atoms, but fails in the molecular analog of 
the Mott regime, when each electron is localized on one 
H atom @,[i]. 

Generalizing to the local spin density approximation 
improves the energy at the cost of an artificial breaking of 
the symmetry. [6| While this can be formally justified, 7] 
in practice it leads to unwanted features. For example 
in the case of a strongly correlated metal the artificial 
breaking of symmetry will give rise to a Fermi surface 
with the wrong Luttinger volume. [8[ 

The breakdown of unpolarized LDA in the H molecule 
can be traced back to the poor treatment of each H atom 
separately. [4j, |5| The density close to the center of the 
H atom is 0.32 a.u. LDA essentially assumes that this 
portion of the system behaves as a uniform electron gas 
with the same density (Wigner-Seitz radius parameter 
r s ~ 0.91). However the two systems have radically dif- 
ferent pair distribution functions which leads to differ- 
ent electron-electron interaction energies (zero for the H 
atom) . LDA thus introduces a spurious interaction of the 
electron with itself, the so-called self-interaction error [9j . 
A way to mitigate this problem would be to have a func- 
tional theory which depends both on the density and the 
pair density (DPDFT) so that it can discriminate be- 
tween situations with the same density but different pair 
densities. 

The idea to involve the pair density in electronic struc- 
ture computations is older than DFT itself lOJjlllI- More 



recently various works explored the possibility to define 
an energy functional based on the pair-density alone [12| - 
Il5| or explore the possibility to adjust the spatially aver- 
aged pair density [16[ but face serious problems on find- 
ing physically acceptable pair densities. [Iff Our proposal 
differs in that we still keep the density as the basic vari- 
able but we use partial information on the pair density 
as an auxiliary variable which gives the functional more 
sensitivity to correlation. In this respect our functional 
is similar to the proposal of Ref. jjj with the difference 
that it does not need an artificial breaking of symmetry. 

To show the feasibility and the usefulness of the formal- 
ism in strongly correlated systems we present a DPDFT 
for the one-dimensional Hubbard model i.e. electrons on 
a lattice with strong local interaction. The functional is 
inspired on the Gutzwiller approximation (GA) [17H19I ] in 
the same way as Kohn-Sham approach is inspired on the 
Hartree approximation. [20j It is thus in principle exact as 
Kohn-Sham theory represents the "exactification" of the 
Hartree approximation [l|. We develop an approximation 
to the unknown functional which becomes exact for uni- 
form systems and involves local or semilocal quantities 
like LDA and its gradient generalizations, but which is 
more accurate and goes beyond LDA in the sense that it 
becomes highly non-local when expressed as a standard 
functional of the density alone. 

We consider a one-dimensional inhomogeneous Hub- 
bard model H — H t + Hjj + H v with 
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here c xa and c xr7 indicate the annihilation and creation 
operator of electrons with spin a on site x and n xa = 
c xa c X(J . t is the nearest-neighbor hopping amplitude 
while U x and v x denote respectively the Hubbard inter- 
action energy and the external potential on site x. For 
reasons to become clear below we allow the interaction 
energy to be site dependent. 

The Hohenberg and Kohn theorem [2l|, |22| guar- 
antees that there exists a functional, ^[{/Ox}] = 
F[{U X , p x }] + J2 X v xPx, which, when minimized with re- 
spect to the density provides the exact ground state en- 



ergy. F[{U X , p x }] is a "universal" functional independent 
of v x , but depending on the U x 's, which represents the 
contribution of H t + Hjj to the energy of a system with 
density p x . A formal DPDFT for this model can be ob- 
tained by performing the Legendre transform, [23| 

T[{d x , Px }} = max ( F[{U X , p x }} - £ U x d x ) (1) 

where d x = (n x ^n x i) is the on-site pair density or double 
occupancy. Identifying the last term in the brackets with 
the interaction energy we arrive at the conclusion that 
T[{d x , p x }] is the interacting kinetic energy of the model 
with the specified pair density and density distributions. 
This is a universal functional which does not depend on 
the specific form of U x nor v x . Below we develop an 
approximation for this functional. 

The ground state energy of the system is obtained by 
minimizing the functional, 



Euv[{d x ,p x }} = T[{d x ,p x }] 
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U x d x + y ' v x p x (2) 



with respect to d x and p x . In the new scheme H t has 
become the fixed part of the Hamiltonian while both Hjj 
and H v are considered as problem dependent. We will 
show below that even in the conventional case in which 
U x is homogeneous the new functional is quite conve- 
nient. 

Eq. © yields, 



dF[{U x ,p x }} 
dU x > 



= d x 



Inversion of this expression determines the set of U x a 
system must have to have the given d x and p x . Since H v 
and Hjj are determined by p x and d x , the wave-function 
and all physical quantities are functionals of p x and d x . 
The interacting kinetic energy can be written as 

T = —t y \Px,x+i + Px+i,x) 

X 

where p x ,x' denotes the one body density matrix, p XtX / = 
^2<j( c x(T C x'a) and is also a functional of p x and d x . 

To proceed we search for a non-interacting system sat- 
isfying the, to be determined, Kohn-Sham like equations, 



/ , i x,x+S l fx+8a ^ 
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and which has the same density of the interacting system. 
Defining the non-interacting one-body density matrix as 
the sum restricted to the occupied Kohn-Sham orbitals, 



P°x,x'= E <pI M W<pP&), 



(4) 
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and the density as p x = p x x , we require thus that, p x 

Px- 



In addition we introduce the kinetic energy reduction 
factors which are also functionals of d and p and satisfy 

Px,x±l = q X ,x±l[{dx> , Px'Y\Px,x±X- ( 5 ) 

With these definitions the kinetic energy functional can 
be written in terms of the non interacting density matrix 
as follows 

<lx>,x>+5[{d x ,pl}]pl, x , +s . (6) 



T[{d x ,p° x }] 



E 
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Similarly to the exchange correlation potential in the 
standard Kohn-Sham approach, q encodes all the com- 
plication of the many-body problem in an exact way. 

Minimizing Eqs. ([2]) with respect to the orbitals and 
using Eq. ([6]) we arrive at 

tx,x±l = t<lx,x±l[{dx> , P X '}], (7) 



..eff 



Eqs. ©-© define the kinetic energy reduction func- 
tional q x , x +s[{d x i , p x '}] and the associated Kohn-Sham 
system. 

For practical computations one needs to introduce 
approximations. In analogy with the Gutzwiller 
approximation 171 llS| we assume the kinetic energy re- 
duction functional factorizes in terms of functions of the 
local densities, q XtV = z(p x , d x )z(p y , d y ). In the following 
we refer to this approximation as "factorized approxi- 
mation" (FA). In the spirit of LDA,[13, 0] we approxi- 
mate the functional dependence of z(p x ,d x ) by request- 
ing that the functional yields the exact energy in the case 
of a uniform system which we obtain in one-dimension 
by solvin g n umerically the exact Bethe-ansatz integral 
equations[25] (in higher dimension a numerical solution 
can be used). Thus z(d,p) is calculated from the con- 
dition z 2 (d,p)T Q (p) = T BA (d,p), where T Q (p) denotes 
the kinetic energy of the non-interacting uniform system. 
T BA (d,p) is obtained as a function of the double occu- 
pancy d : from the exact Bethe-ansatz energy E(U 7 p) by 
the Legendre transform Eq. ([1]) for a system with uni- 
form U and density p. It is also useful to compute the 
function U p (p,d) which yields the Hubbard U a uniform 
system with density p = p x must have, to yield the dou- 
ble occupancy d x of the non-uniform system at the given 
site. We indicate as d(p, U p ) the inverse of U p (p, d). 

The functional Eq. @ has to be minimized with re- 
spect to the local double occupancies d x , and the Kohn- 
Sham orbitals leading to Eq. ((3]). We can use the function 
d(p, U p ) to eliminate d x in favor of a site dependent ef- 
fective interaction, termed the "pseudointeraction" , U p 
with respect to which the minimization is actually done 
[d x — > d(p x ,U p )]. This change of variables allows us to 
avoid the problem of minimizing with respect to a con- 
strained variable. In the case of a homogeneous exter- 
nal potential v x = V, the minimum is attained when 
U p = [/, and one recovers the exact Bethe ansatz result. 




FIG. 1: (a) Error on the energy per site, e, for a periodic 
Hubbard chain of length L = 16 with a binary potential and 
TV = 6 electrons as a function of the potential strength V . 
Notice that as for the LDA, DPDFT-FA becomes asymptot- 
ically exact for a negligible external potential but as the GA 
approximation gives a small error in the case of a strong inho- 
mogeneous potential where LDA fails, (b) Error as a function 
of the interaction strength, (c) Error in the interaction en- 
ergy per site u = U^2 x d x /L as a function of V. (d) The 
pseudointeraction (defined in the text) at odd and even sites 
as a function of the potential strength. 



DPDFT is a variant of DFT. Indeed once a functional 
for T is known we can define a functional of the den- 
sity alone by minimizing over d x , i.e. F[{U, p x }] — 
min/ dx} (T[{d x ,p x }] + UJ2^dx)- Interestingly, as we 
know from previous works[26], this procedure leads to 
a highly non-local functional of the density starting from 
a nearly local functional of variables p x , d x . 

In order to test the functional we solved the DPDFT- 
FA equations for a Hubbard chain with a binary potential 
v x = (— 1) X V (known as ionic Hubbard model[27J, [28|) 
and compare it with standard LDA derived from the 
exact Bethe ansatz solution|24j and with exact results 
obtained by Lanczos diagonalization with the ALPS 
package [29|. 

Fig. Q] shows that DPDFT-FA performs much better 
than LDA when the system becomes highly inhomo- 
geneous (a) and strongly interacting (b). The double 
occupancy in LDA is independent of the environment 
which leads to a poor approximate interaction energy for 
strongly inhomogeneous systems (c). In DPDFT-FA the 
double occupancy is allowed to adapt, so that localized 
electrons, for large values of the potential, tend to have 
a small double occupancy and a small interaction energy 
(c), explaining the better performance of DPDFT-FA re- 
spect to LDA. The enhanced pseudointeraction on the 
more charged sites [odd x in (d)] leads to a reduced dou- 
ble occupancy and explains the small interaction energy. 

The behavior of DPDFT-FA is qualitatively similar to 



the GA except at small V where the GA is obviously 
not exact. The errors in other quantities like kinetic en- 
ergy, potential energy and density (not shown) are also 
generically smaller in DPDFT-FA than in LDA. 

These results suggest that the self-interaction (SI) 
should be strongly reduced in DPDFT-FA with respect 
to LDA. In order to verify that this is the case we have 
solved the problem of one electron with one attractive 
impurity site in a chain. The potential is given by 
v x = — V5 x fi. This is the lattice analogue of the hy- 
drogen atom in the continuum. Being a single electron 
problem, the exact solution has d x ~ for all U while 
both LDA and DPDFT-FA yield a finite d x . We define 
the self-interaction error as the spurious interaction en- 
ergy of the single electron problem, Esi — U J2 X d x ■ 

In Fig. Ufa) we plot the self-interaction error as a func- 
tion of the potential strength at the impurity site, with 
interaction parameter U = 4t. For large V the charge be- 
comes localized at the impurity site with po ~ 1 • I n LDA 
the interaction energy corresponds to that of a nearly 
half-filled uniform Hubbard model which is clearly a very 
bad approximation thus E$i (red dashed line) increases 
and tends to saturate at a large value. In DPDFT- 
FA (blue solid line) E$i starts with a slower increase 
and then decays to very small values. In this case the 
total double occupancy Esi/U becomes small showing 
the adaptability of the pair density to the local envi- 
ronment. As shown in Fig. H^b), the reduction of the 
self-interaction error is large for a wide range of the in- 
teraction. 

It is interesting to compare DPDFT with traditional 
Kohn-Sham DFT for the same model 0. In the lat- 
ter the difference between the interacting and the non- 
interacting kinetic energy is absorbed additively in the 
exchange-correlation potential. Here the correction is 
multiplicative and included in the kinetic energy reduc- 
tion factor functional q. 




FIG. 2: Self-interaction error for a chain of L — 12 sites, a sin- 
gle electron and a potential in which one site has strength —V 
as a function of V (a) and U (b). For large potential strength 
the density at the impurity becomes close to one. In LDA 
this fixes the double occupancy to an unphysical value while 
in DPDFT-FA the double occupancy adjusts to small values 
when the electron localizes leading to a small self-interaction. 



Formally the minimization in DFT has to be restricted 
to densities, and here also double occupancies, that cor- 
respond to a physical wave function (N-representability). 
While this may appear as a severe difficulty [15] it is of- 
ten not a problem in practical implementations. It has 
indeed not hampered the development of DFT methods. 
In our case the problem is not more severe than in LDA 
because each portion of the system is approximated by a 
uniform system but with a modified interaction. 

We have shown DPDFT at work in the lattice but a 
similar functional can be defined in the continuum. Ide- 
ally one would like to use the pair density, 7 CT(T '(r, r') = 

("0i( r )Vv ( r ')Vv ( r ')V ; (T( r )) ) as a variable, with Vv(r) the 
field operator at point r with spin a. This, however, 
would be quite cumbersome as in each point a full func- 
tion must be determined. A more practical approach is to 
impose a spatial dependent constraint on the pair den- 
sity and use such a constraint to parameterize families 
of physical pair density functions. For example we can 
define, 



D(v) 



E 



d 3 r' 1<T(T ,(r,r')d(a-\r-r'\), 



with the Heaviside function and a is an appropri- 
ately chosen cutoff radius. D measures the double oc- 
cupancy probability within a sphere of radius a. These 
or other constraints can be implemented 23] by replacing 
the physical interaction w(r,r') with a fictitious inter- 
action. In the present example the fictitious interaction 
would read w(r, r') + U(r)9(a — |r — r'Q. As in the lat- 
tice, the corresponding Hohenberg and Kohn functional 
has U(r) as a variable, F[n(r), U(r)] which allows to de- 
fine the functional 



G[n(r),D(r)} 



max 

U(r) 



F[n{r),U(r)} 



J d 3 rU(r)D(r] 



leading to a theory where both n(r) and D(r) are fun- 
damental variables similar to the lattice but with the 
difference that G is not the kinetic energy. Solution of 
the uniform problem in the presence of the fictitious in- 
teraction with a constant U(r) = U p could serve as a 
basis for approximate functionals which converge to the 
LDA in the uniform case but have an adaptive exchange 
correlation hole in non-uniform situations. 

We have conceptually shown how a DFT which uses 
the pair density as an auxiliary variable can be intro- 
duced and we have developed an approximation for the 
Hubbard model inspired on the GA combined with LDA 
ideas. Formally DPDFT can also be defined in the con- 
tinuum with a local variable D(r) playing the role of the 
double occupancy in the lattice. Approximate function- 
als based on this approach should allow more control on 
the correlations built on the underlying wave- function re- 
spect to what LDA does, and could help to extend the 



success of DFT methods to strongly correlated systems 
where correlations are substantially different from those 
of the homogeneous electron gas. 
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